package scu.maqiang.numeric;

import scu.maqiang.csparsej.tdouble.*;

import java.util.Arrays;

/**
 * 稀疏矩阵直接求解器，包装调用了csparsej包中的直接求解器, 可实现实系数与复系数稀疏矩阵直接求解
 * @author 马强
 *
 */
public class DirectSSolver extends Solver{
	Dcs_common.Dcs matrix = null;
	Dcs_common.Dcss S = null;
	Dcs_common.Dcsn N = null;
	int n = 0;

	public DirectSSolver() {
	}

	/**
	 * 构造函数, 指定求解器要操作的实稀疏矩阵
	 * @param mat 实稀疏矩阵
	 */
	public DirectSSolver(SRMatrix mat) {
		Dcs_common.Dcs temp = mat.toDcs_matrix();
		matrix = Dcs_compress.cs_compress(temp);
		n = mat.getM();
	}


	/**
	 * 构造函数, 指定求解器要操作复稀疏矩阵
	 * @param mat 复稀疏矩阵
	 */
	public DirectSSolver(SCMatrix mat) {
		int n = mat.getN();
		this.n = 2 * n;
		SRMatrix ABMatrix = new SRMatrix(n * 2);
		for(int i = 0; i < n; i++) {
			for(SCEntry ele: mat.mat[i]) {
				ABMatrix.setElementFast(i, ele.col, ele.val.real);
				ABMatrix.setElementFast(i, ele.col + n, -ele.val.image);
				ABMatrix.setElementFast(i + n, ele.col, ele.val.image);
				ABMatrix.setElementFast(i + n, ele.col + n, ele.val.real);
			}
		}
		Dcs_common.Dcs temp = ABMatrix.toDcs_matrix();
		matrix = Dcs_compress.cs_compress(temp);
	}

	public void setMatrix(SRMatrix mat) {
		Dcs_common.Dcs temp = mat.toDcs_matrix();
		matrix = Dcs_compress.cs_compress(temp);
		n = mat.getM();
	}

	public void setMatrix(SCMatrix mat) {
		int n = mat.getN();
		this.n = 2 * mat.getM();
		SRMatrix ABMatrix = new SRMatrix(n * 2);
		for(int i = 0; i < n; i++) {
			for(SCEntry ele: mat.mat[i]) {
				ABMatrix.setElementFast(i, ele.col, ele.val.real);
				ABMatrix.setElementFast(i, ele.col + n, -ele.val.image);
				ABMatrix.setElementFast(i + n, ele.col, ele.val.image);
				ABMatrix.setElementFast(i + n, ele.col + n, ele.val.real);
			}
		}
		Dcs_common.Dcs temp = ABMatrix.toDcs_matrix();
		matrix = Dcs_compress.cs_compress(temp);
	}

	/**
	 * 使用LU分解方法求解一般系数矩阵的线性方程组
	 * @param RHS 线性方程组右端项
	 * @param order 使用的排序方法 (0 to 3)
	 * @param tol 部分选主元的容许误差
	 * @return 方程组的解
	 */
	public double[] LUSolve(double[] RHS, int order, double tol) {
		double[] result = RHS.clone();
		long begin = System.nanoTime();
		Dcs_lusol.cs_lusol(order, matrix, result, tol);
		long end = System.nanoTime();
		System.out.printf("min: %20.10e\t max: %20.10e\n", MVO.min(result), MVO.max(result));
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		return result;
	}
	/**
	 * 使用Cholesky分解方法求解对称正定矩阵的线性方程组
	 * @param RHS 线性方程组右端项
	 * @param order 使用的排序方法 (0 to 3)
	 * @return 方程组的解
	 */
	public double[] CholeskySolve(double[] RHS, int order) {
		double[] result = RHS.clone();
		long begin = System.nanoTime();
		Dcs_cholsol.cs_cholsol(order, matrix, result);
		long end = System.nanoTime();
		System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(result), MVO.max(result));
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		return result;
	}

	public void CholeskySolve(double[] RHS, double[] result, int order) {
		double[] tempResult = RHS.clone();
		ParamCheck.checkEqual(result.length, RHS.length);
		long begin = System.nanoTime();
		Dcs_cholsol.cs_cholsol(order, matrix, tempResult);
		long end = System.nanoTime();
		System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(result), MVO.max(result));
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		Arrays.setAll(result, i -> tempResult[i]);
	}


	/**
	 * QR方解方法求解最小二乘问题
	 * @param RHS 最小二乘右端项
	 * @param order 使用的排序方法 (0 to 3)
	 * @return 线性方程组最小二乘解
	 */
	public double[] QRSolve(double[] RHS, int order) {
		double[] result = RHS.clone();
		Dcs_qrsol.cs_qrsol(order, matrix, result);
		return result;
	}

	/**
	 * 使用LU分解方法求解复系数线性方程组
	 * @param RHS 复线性方程组右端项
	 * @param order 使用的排序方法 (0 to 3)
	 * @param tol 选主元容许误差
	 * @return 复数解数组
	 */
	public Complex[] complexSolve(Complex[] RHS, int order, double tol) {
		int n = RHS.length;
		ParamCheck.checkEqual(matrix.m, 2 * n);
		double[] tempResult = new double[2 * n];
		for(int i = 0; i < n; i++) {
			tempResult[i] = RHS[i].real;
			tempResult[i + n] = RHS[i].image;
		}
		long begin = System.nanoTime();
		Dcs_lusol.cs_lusol(order, matrix, tempResult, tol);
		System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(tempResult), MVO.max(tempResult));
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		Complex[] result = new Complex[n];
		for(int i = 0; i < n; i++) {
			result[i] = new Complex(tempResult[i], tempResult[i + n]);
		}
		return result;
	}

	/**
	 * 对对称正定矩阵进行Cholesky分解, 与Solve_Cholesky一起使用
	 * @param order 排序方法 (0 to 3)
	 * @see #Solve_Cholesky(double[], int)
	 * @return 是否分解成功
	 */
	public boolean Factorize_Cholesky(int order) {
		boolean ok;
		n = matrix.n;
		if (!Dcs_util.CS_CSC(matrix))
			return (false); /* check inputs */
		S = Dcs_schol.cs_schol(order, matrix); /* ordering and symbolic analysis */
		N = Dcs_chol.cs_chol(matrix, S); /* numeric Cholesky factorization */
		ok = (S != null && N != null);
		//System.out.println(Arrays.toString(b));
		return (ok);
	}

	/**
	 * 对已经分解好的对称正定矩阵进行求解, 与Factorize_Cholesky一起使用
	 * @param b 线性方程组右端项
	 * @see #Factorize_Cholesky(int)
	 * @return 方程组解
	 */
	public double[] Solve_Cholesky(double[] b, int verbosity) {
		double x[] = new double[n];
		double bb[] = b.clone();
		boolean ok = (S != null && N != null && b != null);
		//System.out.println(ok);
		long begin = System.nanoTime();
		if (ok) {
			Dcs_ipvec.cs_ipvec(S.pinv, bb, x, n); /* x = P*b */
			Dcs_lsolve.cs_lsolve(N.L, x); /* x = L\x */
			Dcs_ltsolve.cs_ltsolve(N.L, x); /* x = L'\x */
			Dcs_pvec.cs_pvec(S.pinv, x, bb, n); /* b = P'*x */
		}
		long end = System.nanoTime();
		if(verbosity >= 1) {
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(bb), MVO.max(bb));
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		}
		return bb;
		//System.out.println(Arrays.toString(b));
	}

	public void Solve_Cholesky(double[] b, double[] xx, int verbosity) {
		double x[] = new double[n];
		double bb[] = b.clone();
		boolean ok = (S != null && N != null && b != null);
		ParamCheck.checkEqual(xx.length, n);
		//System.out.println(ok);
		long begin = System.nanoTime();
		if (ok) {
			Dcs_ipvec.cs_ipvec(S.pinv, bb, x, n); /* x = P*b */
			Dcs_lsolve.cs_lsolve(N.L, x); /* x = L\x */
			Dcs_ltsolve.cs_ltsolve(N.L, x); /* x = L'\x */
			Dcs_pvec.cs_pvec(S.pinv, x, bb, n); /* b = P'*x */
		}
		long end = System.nanoTime();
		if(verbosity >= 1) {
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(bb), MVO.max(bb));
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		}
		Arrays.setAll(xx, i-> bb[i]);
		//System.out.println(Arrays.toString(b));
	}

	/**
	 * 对一般系数矩阵进行LU分解, 与Solve_LU一起使用
	 * @param order 排序方法 (0 to 3)
	 * @see #Solve_LU(double[])
	 * @return 是否分解成功
	 */
	public boolean Factorize_LU(int order, double tol) {
		int n;
		boolean ok;
		if (!Dcs_util.CS_CSC(matrix))
			return (false); /* check inputs */
		S = Dcs_sqr.cs_sqr(order, matrix, false); /* ordering and symbolic analysis */
		N = Dcs_lu.cs_lu(matrix, S, tol); /* numeric LU factorization */
		ok = (S != null && N != null);
		//System.out.println(Arrays.toString(b));
		return (ok);
	}

	/**
	 * 对已经LU分解好的系数矩阵进行求解, 与Factorize_LU一起使用
	 * @param b 线性方程组右端项
	 * @see #Factorize_LU(int, double)
	 * @return 方程组解
	 */
	public double[] Solve_LU(double[] b, int verbosity) {
		int n = matrix.n;
		double[] x = new double[n]; /* get workspace */
		double bb[] = b.clone();
		boolean ok = (S != null && N != null && b != null);
		long begin = System.nanoTime();
		if (ok) {
			Dcs_ipvec.cs_ipvec(N.pinv, bb, x, n); /* x = b(p) */
			Dcs_lsolve.cs_lsolve(N.L, x); /* x = L\x */
			Dcs_usolve.cs_usolve(N.U, x); /* x = U\x */
			Dcs_ipvec.cs_ipvec(S.q, x, bb, n); /* b(q) = x */
		}
		long end = System.nanoTime();
		if(verbosity >= 1) {
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(bb), MVO.max(bb));
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		}
		return bb;
	}

//	/**
//	 * 当使用这个构造函数时，求解问题使用PPLUSolve
//	 * @param A
//	 */
//	public DirectSSolver(SRMatrix A) {
//		this.A = A;
//		int m = A.getM();
//		P = A.symamd();
//		L = new SRMatrix(m);
//		U = new SRMatrix(m);
//		//P.setDiag(1.0);
//		L.setDiag(1.0);
//		LU();
//	}
//
//	/**
//	 *  当使用这个构造函数时，求解问题使用PLUSolve
//	 * @param A
//	 */
//	public DirectSSolver(SRMatrix A, int nBand) {
//		this.A = A;
//		this.nBand = nBand;
//		int m = A.getM();
//		//P = new SRMatrix(m);
//		L = new SRMatrix(m);
//		U = new SRMatrix(m);
//		//P.setDiag(1.0);
//		L.setDiag(1.0);
//		P = new int[m];
//		Arrays.setAll(P, i -> i);
//		PLU();
//		//System.out.println(P);
//	}
//
//	/**
//	 * 对稀疏矩阵A进行LU分解
//	 */
//	public void LU() {
//		long begin = System.nanoTime();
//		int m = A.getM();
//		int fillinCount = 0;
//		for(int i = 0; i < m - 1; i++) {
//			double maxval = A.getElement(i, i);
//			for(int j = i + 1; j < m; j++) {
//				SREntry ji = A.getSREntry(j, i);
//				if (ji != null) {
//					double ratio = ji.val / maxval;
//					ji.val = ratio;
//					for(SREntry ele : A.mat[i]) {
//						if (ele.col > i) {
//							SREntry jCol = A.getSREntry(j, ele.col);
//							if(jCol != null) {
//								jCol.val -= ele.val * ratio;
//							} else {
//								A.setElementFast(j, ele.col, -ele.val * ratio);
//								System.out.println("Fill in in\t " + j + "\t, " + ele.col);
//								fillinCount++;
//							}
//						}
//					}
//				}
//			}
//		}
//		System.out.println("Factorizing finished....");
//		System.out.println("Fill in number: " + fillinCount);
//		long end = System.nanoTime();
//		System.out.println("Factorizing elapsed...\t" + (end - begin) / 1.0e9 + " s");
//		//System.out.println("P:" + Arrays.toString(P));
//		//System.out.println("P:" + MVO.toString(P.toDenseMatrix()));
////		System.out.println("A:" + MVO.toString(A.toDenseMatrix()));
//		for(int i = 0; i < m; i++) {
//			for(SREntry ele : A.mat[i]) {
//				if (ele.col >= i) {
//					U.mat[i].add(new SREntry(ele.col, ele.val));
//				} else {
//					L.mat[i].add(new SREntry(ele.col, ele.val));
//				}
//			}
//		}
//	}
//
//	public void PPLUSolve(double[] RHS, double[] x, int verbosity) {
//		long begin = System.nanoTime();
//		int m = L.getM();
//		double[] Pb = new double[m];
//		Arrays.setAll(Pb, i -> RHS[P[i]]);
//		double[] C = new double[m];
//		C[0] = Pb[0];
//		for(int i = 1; i < m; i++) {
//			double temp = 0.0;
//			for(SREntry ele : L.mat[i]) {
//				if (ele.col < i) {
//					temp += C[ele.col] * ele.val;
//				}
//			}
//			C[i] = Pb[i] - temp;
//		}
//
//		double[] xx = new double[m];
//		double vMM = U.getElement(m - 1, m - 1);
//		xx[m - 1] = C[m - 1] / vMM;
//		for(int i = m - 2; i >=0; i--) {
//			double temp = 0.0;
//			for(SREntry ele : U.mat[i]) {
//				if (ele.col > i) {
//					temp += xx[ele.col] * ele.val;
//				} else {
//					vMM = ele.val;
//				}
//			}
//			xx[i] = (C[i] - temp) / vMM;
//		}
//		for(int i = 0; i < m; i++) {
//			x[P[i]] = xx[i];
//		}
//
//		if (verbosity >= 1) {
//			System.out.printf("min: %20.10e\t max: %20.10e\n", MVO.min(x), MVO.max(x));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//	}
//
//
//
//	/**
//	 * 对稀疏矩阵A进行PLU分解
//	 */
//	public void PLU() {
//		long begin = System.nanoTime();
//		int m = A.getM();
//		int fillinCount = 0;
//		for(int i = 0; i < m - 1; i++) {
//			Pair<Integer, Double> maxIdxResult = A.maxAbsColValDown(i, nBand);
//			int idxJ = maxIdxResult.getFirst();
//			double maxval = maxIdxResult.getSecond();
//			if (Math.abs(maxval) < Constants.Er) {
//				throw new ArithmeticException("Zero pivot element!");
//			}
//			if (maxIdxResult.getFirst() != i) {
//				A.exchangeRow(i, idxJ);
//				//System.out.println("OK!");
//				//P.exchangeRow(i, idxJ);
//				int temp = P[i];
//				P[i] = P[idxJ];
//				P[idxJ] = temp;
//			}
//			int end = Math.min(i + nBand + 3, m);
//
//			for(int j = i + 1; j < end; j++) {
//				SREntry ji = A.getSREntry(j, i);
//				if (ji != null) {
//					double ratio = ji.val / maxval;
//					ji.val = ratio;
//					for(SREntry ele : A.mat[i]) {
//						if (ele.col > i) {
//							SREntry jCol = A.getSREntry(j, ele.col);
//							if(jCol != null) {
//								jCol.val -= ele.val * ratio;
//							} else {
//								A.setElementFast(j, ele.col, -ele.val * ratio);
//								fillinCount++;
//								System.out.println("Fill in in\t " + j + "\t, " + ele.col);
//							}
//						}
//					}
//				}
//			}
//		}
//		System.out.println("Factorizing finished....");
//		System.out.println("Fill in number: " + fillinCount);
//		long end = System.nanoTime();
//		System.out.println("Factorizing elapsed...\t" + (end - begin) / 1.0e9 + " s");
//		//System.out.println("P:" + Arrays.toString(P));
//		//System.out.println("P:" + MVO.toString(P.toDenseMatrix()));
////		System.out.println("A:" + MVO.toString(A.toDenseMatrix()));
//		for(int i = 0; i < m; i++) {
//			for(SREntry ele : A.mat[i]) {
//				if (ele.col >= i) {
//					U.mat[i].add(new SREntry(ele.col, ele.val));
//				} else {
//					L.mat[i].add(new SREntry(ele.col, ele.val));
//				}
//			}
//		}
//
////		System.out.println(Arrays.toString(P));
////		System.out.println("P:" + MVO.toString(P.toDenseMatrix()));
////		System.out.println("L:" + MVO.toString(L.toDenseMatrix()));
////		System.out.println("U:" + MVO.toString(U.toDenseMatrix()));
//	}
//
//	/**
//	 *
//	 * @param RHS
//	 * @param x
//	 * @param verbosity
//	 */
//	public void PLUSolve(double[] RHS, double[] x, int verbosity) {
//		long begin = System.nanoTime();
//		int m = L.getM();
//		double[] Pb = new double[m];
//		Arrays.setAll(Pb, i -> RHS[P[i]]);
//		double[] C = new double[m];
//		C[0] = Pb[0];
//		for(int i = 1; i < m; i++) {
//			double temp = 0.0;
//			for(SREntry ele : L.mat[i]) {
//				if (ele.col < i) {
//					temp += C[ele.col] * ele.val;
//				}
//			}
//			C[i] = Pb[i] - temp;
//		}
//
//		double vMM = U.getElement(m - 1, m - 1);
//		x[m - 1] = C[m - 1] / vMM;
//		for(int i = m - 2; i >=0; i--) {
//			double temp = 0.0;
//			for(SREntry ele : U.mat[i]) {
//				if (ele.col > i) {
//					temp += x[ele.col] * ele.val;
//				} else {
//					vMM = ele.val;
//				}
//			}
//			x[i] = (C[i] - temp) / vMM;
//		}
//
//		if (verbosity >= 1) {
//			System.out.printf("min: %20.10e\t max: %20.10e\n", MVO.min(x), MVO.max(x));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//	}
//
//
//	SRMatrix A;
//	int[] P;
//	SRMatrix L;
//	SRMatrix U;
//	int nBand;
}
